#include "PLinearMotion.h"

PLinearMotion::PLinearMotion(const PLinearMotion& p0)
{
	dragForce = vector<double>(2, 0.0); accDrag = vector<double>(2, 0.0);

	Rep = p0.Rep, coefDrag = p0.coefDrag;
	dragForce[0] = p0.dragForce[0], dragForce[1] = p0.dragForce[1];
	accDrag[0] = p0.accDrag[0], accDrag[1] = p0.accDrag[1];
}

PLinearMotion::PLinearMotion(const Slip& pSlip, const Particle& p_in, const Cphase& cphase_in)
{
	dragForce = vector<double>(2, 0.0); accDrag = vector<double>(2, 0.0);

	Rep = pSlip.Rep;
	coefDrag = Cd(Rep);
	double tao = 4.0 / 3.0 * p_in.rho * p_in.diam * p_in.diam / (Rep * coefDrag * cphase_in.mu);
	dragForce[0] = M_PI / 8.0 * coefDrag * p_in.diam * p_in.diam * 
		abs(pSlip.slipVelocity[0]) * pSlip.slipVelocity[0] * cphase_in.rho;
	dragForce[1] = M_PI / 8.0 * coefDrag * p_in.diam * p_in.diam *
		abs(pSlip.slipVelocity[1]) * pSlip.slipVelocity[1] * cphase_in.rho;

	accDrag[0] = 0.75 * coefDrag * abs(pSlip.slipVelocity[0]) * pSlip.slipVelocity[0] * cphase_in.rho
		/ (p_in.diam * p_in.rho);
	accDrag[1] = 0.75 * coefDrag * abs(pSlip.slipVelocity[1]) * pSlip.slipVelocity[1] * cphase_in.rho
		/ (p_in.diam * p_in.rho);

}
double PLinearMotion::Cd(double Rep)
{
	// empirical spherical model
	double Cd = 0.0;
	if (Rep < 0.1)
	{
		Cd = 24.0 / Rep;
	}
	else if (Rep >= 0.1 && Rep < 1.0)
	{
		Cd = 22.73 / Rep + 0.0903 / (Rep * Rep) + 3.69;
	}
	else if (Rep >= 1.0 && Rep < 10.0)
	{
		Cd = 29.1667 / Rep - 3.8889 / (Rep * Rep) + 1.222;
	}
	else if (Rep >= 10.0 && Rep < 100.0)
	{
		Cd = 46.5 / Rep - 116.67 / (Rep * Rep) + 0.6167;
	}
	else if (Rep >= 100.0 && Rep < 1000.0)
	{
		Cd = 98.33 / Rep - 2778.0 / (Rep * Rep) + 0.3644;
	}
	else if (Rep >= 1000.0 && Rep < 5000.0)
	{
		Cd = 148.62 / Rep - (4.75e4) / (Rep * Rep) + 0.357;
	}
	else if (Rep >= 5000.0 && Rep < 10000.0)
	{
		Cd = -490.546 / Rep + (57.87e4) / (Rep * Rep) + 0.46;
	}
	else if (Rep >= 10000.0 && Rep < 50000.0)
	{
		Cd = -1662.5 / Rep + (5.4167e6) / (Rep * Rep) + 0.5191;
	}
	else
	{
		assert(false && "Cd value is out of Range!");
	}
	return Cd;
}